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Thermodynamics of the t-J Ladder: A Stable Finite Temperature Density Matrix 

Renormahzation Group Calculation 
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Accurate numerical simulations of a doped t-J model on a two-leg ladder are presented for the par- 
ticle number, chemical potential, magnetic susceptibility and entropy in the limit of large exchange 
coupling on the rung using a finite temperature density matrix renormahzation group (TDMRG) 
method. This required an improved algorithm to achieve numerical stability down to low tempera- 
tures. The thermal dissociation of hole pairs and of the rung singlets are separately observed and 
the evolution of the hole pair binding energy and magnon spin gap with hole doping is determined. 
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Standard quantum Monte Carlo methods for the sim- 
ulation of fcrmions are limited to relatively high tem- 
peratures due to the fermion sign problem. The density 
matrix renormahzation group method (DMRG) [Q allows 
simulations of large clusters but is limited to groundstate 
properties. In this Letter we report on a finite tempera- 
ture DMRG (TDMRG) f||| method using an improved 
numerically stable algorithm to simulate a strongly in- 
teracting fermion system down to low temperatures. 

The TDMRG method applies the DMRG to the quan- 
tum transfer matrix (QTM) in the real space direction 
[§. In the TDMRG iterations the QTM is enlarged in 
the imaginary time direction and iterates to lower tem- 
peratures at fixed Trotter time steps At. This is in con- 
trast to the DMRG method in which the system grows in 
the real space direction. The TDMRG has the advantage 
that the free energy and other thermodynamic quantities 
for the infinite system can be obtained directly from the 
largest eigenvalue and the QTM and of the corresponding 
eigenvector. 

The system we examine is a two-leg t-J ladder model in 
the limit where the exchange interaction across the rungs 
(J 7 ) is large compared to the value along the legs( J) and 
to the isotropic hopping integral t. The ground state 
properties of this model at low hole doping have been an- 
alyzed previously by exact diagonalization of small clus- 
ters . In this limit J' 3> </, t the thermal dissociation 
of hole pairs and the excitation of triplet magnons can 
be distinguished. These strong coupling processes are a 
good test for any method. 

In this Letter we present accurate results for the mag- 
netic susceptibility \i the particle number n and the en- 
tropy density s in the grand canonical ensemble as a func- 
tion of chemical potential /i and temperature T and then 
remap x(/i, T) —* x( n : T) to obtain the T-dependence at 
constant density. 

Previous versions of the TDMRG method for fermions 
H have suffered from numerical instabilities due to the 
non-Hcrmiticity of the QTM and the corresponding den- 
sity matrices which are constructed from the right and 



left eigenvector of the largest eigenvalue of the QTM. 
These numerical instabilities grow as the number of 
states kept is increased or the filling is changed away 
from half-filled bands. We have identified the loss of 
biorthonormality between the left and right eigenvectors 

(Vi,Vj) — Sij of the density matrix as the source of 
the problem. 
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become biorthonormal, which leads to severe loss of pre- 
cision due to roundoff errors if the overlap between these 
vectors is small. These near-breakdowns occur especially 
often in conjunction with the second numerical problem, 
spurious small imaginary parts of (nearly) degenerate 
eigenvalue pairs. This latter problem can be solved by us- 
ing the real and imaginary components of the correspond- 
ing complex conjugate eigenvector pairs and discarding 
the imaginary part of the eigenvalues, which are artifacts 
of roundoff errors and are only of the order of the ma- 
chine precision. To circumvent the loss of precision in the 
former problem of nearly orthogonal eigenvectors our al- 
gorithm uses an iterative re-biorthogonalization step (7j 
for the eigenvectors kept, which stabilizes the method 
for all temperatures. Technical details of the algorithms 
will be presented elsewhere ||. In Tab. | we show re- 
sults of numerical stability tests of the original and 
our improved algorithm for the case of noninteracting 
spin 5=1/2 fermions in one dimension. For this simple 
fermionic model the original TDMRG method becomes 
numerically unstable whenever more than about m = 10 
states are kept, thus severely restricting the achievable 
accuracy. The improved algorithm presented here, on 
the other hand, is always numerically stable and achieves 
much higher accuracy. The test example clearly demon- 
strates the need for numerical stabilization in the simu- 
lation of fermionic models. 

The results of the stabilized TDMRG method are accu- 
rate and unbiased, with errors only originating from the 
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TABLE I. Free energy density / for noninteracting spin 
S — 1/2 fermions on a chain with /j = at a temperature 
T = 0.1 after a hundred DMRG steps (At = 0.1), where 
energies are given in units of the hopping integral t. The 
first column (I) is the original TDMRG algorithm ||[|, and 
the second column (II) is our improved method including the 
re-biorthogonalization step. Cases where the algorithm be- 
comes numerically unstable are denoted with tn, where n is 
the number of DMRG steps that could be performed success- 
fully. 
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finite size of the Trotter time steps and the truncation in 
the DMRG algorithm. The latter are usually very small 
if the number of states kept, to, is large enough, and the 
former can be eliminated by extrapolating At — > by fit- 
ting to a polynomial in At 2 . We have used Trotter time 
steps from Art = 0.01 to Art = 0.2, and to between 
to = 40 and to = 60. We made use of spin conservation 
symmetry, the subspace of zero winding number and the 
reflection symmetry of the ladder along the rungs to op- 
timize the calculations and reduce numerical errors. 

Thermodynamic quantities such as the internal energy 
U, the hole density n/t(= 1 — n) and the magnetic suscep- 
tibility x have been determined directly from the eigen- 
vectors of the transfer matrix Q|| . This is preferable to 
taking numerical derivatives of the free energy density 
obtained from the largest eigenvalue of the QTM. 

The low temperature properties of a doped t-J two-leg 
ladder in the limit J' 3> J, * are determined solely by the 
singlet hole pairs (HP) ||. They form a hard core boson 
gas with a bandwidth of At*, with t* = 2* 2 ( J' -At 2 /J')' 1 
in second order perturbation theory. Neglecting a weak 
nearest neighbor attraction, the HP fluid can be mapped 
to an ideal Fermi gas in this one dimensional geometry. 
As the temperature T is increased the HPs dissociate into 
two quasiparticles (QP), each consisting of a single elec- 
tron with spin S = 1/2 in a rung bonding state. Each 
QP propagates with a bandwidth of It so that in the 
limit of low hole doping n./, <C 1 the HP binding energy 
is E B = J 1 - At + 4* 2 (J' - At 2 / J')- 1 . The gas of QPs 
with density uqp(T) contributes to the spin susceptibil- 
ity as a nondegenerate gas of S = 1/2 fermions. A second 
contribution comes from the thermal excitation of singlet 
rungs to a triplet magnon state. The activation energy 
for a magnon Am (= J' — J+J 2 /2J' in second order per- 
turbation theory in the limit n/, — » [[l0|]) is larger than 
that of the QPs (Aqp = Eb/2) but since the density of 



FIG. 1. Hole density nn as a function of electron chemical 
potential fj, and temperature T in the strong coupling regime 
J = t/2 — J'/IO. The dashed lines are fits to a hard core bo- 
son model for the hole pairs. Note, the HP chemical potential 
is — 2/i. 

QPs is limited by the hole density (riQp(T) < rih), the 
temperature evolution of x( n > T) is determined largely 
by the magnons at low doping. 

We now turn to the presentation of our finite tem- 
perature results obtained using the improved TDMRG 
algorithm for J = t/2 = J'/IO and compare them to ex- 
pectations based on the above discussion of this strong 
coupling regime. As the calculations were performed in 
the grand canonical ensemble we first present results for 
the hole density n.fc(/x, T). A selection of our results is 
presented in Fig. ^ including a fit to a hard core boson 
model for the HPs: 

e k H p = e H p + 2t*cosk + 2fj,n h . (1) 

Fitting the data for < 0.1 at temperatures T < 0.5* 
we obtain an estimate for the center of the band for HPs 
at e H p = 4.82(6)* and a bandwidth At* = 1.5(2)*. The 
minimum energy to add a HP to an undoped ladder is 
€rp — 2t* = 4.1(1)*, in good agreement with values from 
the finite clusters (e H p ~ 4.71*, At* w 1.494* @). 

A further confirmation of the validity of this hard core 
boson model for the HPs comes from considering the low 
temperature entropy density per site s, determined from 
the free energy density / and the energy density u as 
s — (u — f)/T. As can be seen in Fig. |] the entropy at 
T < 0.3* and low doping (n.^ < 0.1) is also well described 
by the hard core boson model for the HPs. 

At higher temperatures the thermal dissociation of 
HPs into two independent QPs and the thermal exci- 
tation of magnons from rung singlets govern the thermo- 
dynamics. These processes show up in the spin suscep- 
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FIG. 2. Entropy density s as a function of chemical po- 
tential n and temperature T in the strong coupling regime 
J = t/2 — J' / 10. The dashed lines are the values for the 
same hard core boson model as in Fig. [l], using the parame- 
ters obtained in that fit. 



FIG. 3. Uniform magnetic susceptibility per site %t °f the 
t-J ladder for J = t/2 — J'/IO and different hole-densities n/,. 
The symbols denote the results of the TDMRG algorithm, and 
the solid lines are the fitted curves according to Eq. (|^) . The 
fitting parameters are listed in Tab. |n|. 



tibility x(T), which is easiest to interpret in the canon- 
ical ensemble with fixed hole density nh- Therefore we 
use rih(n,T) to remap x(PiT) — > x( n h,T). The val- 
ues of x(a*j T) were calculated by measuring the magne- 
tization (S Z (T)) in the presence of a small external field 
h/t = 5 x 10~ 3 . The results for x( n h, T) appear in Fig. [|. 

At high temperatures T J' , x follows a Curie- 
law for free spins X = 0- ~ n>h)/4T, and it decreases 
when the temperature is lowered below the magnon-gap 
Am ~ 4.13t. The maximum of the peak is shifted to- 
wards lower T with increasing doping, indicating a re- 
duction of the magnon gap due to interactions with holes. 
Simultaneously the magnon bandwidth is enhanced, in- 
dicating that the energy of a localized magnon is not 
much changed by the holes. At very low temperatures of 
T < 0.5t we can see a second exponential decrease of x 
with a smaller gap, which we attribute to the recombina- 
tion of QPs into HPs at temperatures below the QP-gap, 
Aqp = Eb/2. Note the magnitude of this contribution 
increases with rih. 

A quantitative description of x( n h, T) can be given by 
adding separately the contributions of the QPs xqp an d 
of the magnons xm, i.e.: 



x(nh,T) = XQp(nh,T) + x,M(n h ,T). 



(2) 



Xqp is approximated by the value for free spins Xqp — 
uqp{T) / AT with a temperature dependent density of the 
QPs determined by the energy dispersion of the QPs 
Aqp + a QP {l + cos fc)/2 with = 1/T: 
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n QP = v / dk ~0^ k 



(3) 



TABLE II. Gap of the spin S = 1/2 quasi-particles 
Aqp and magnon gap Am, as well as the parameters oqp 
(<2m) which determine the bandwidth of the quasi-particles 
(magnons) obtained by fitting Eq. (^|) to our TDMRG data 
for different hole densities rih- 
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Aqp 


A M 


aQP 


a m 


0.0 




4.1(1) 




0.7(1) 


0.025 


0.7(1) 


3.4(1) 


0.9(2) 


1.6(2) 


0.05 


0.8(1) 


3.3(1) 


0.9(2) 


1.8(2) 


0.1 


1.0(1) 


3.3(1) 


0.6(3) 


1.7(2) 


0.15 


0.9(1) 


3.2(1) 


1.3(2) 


2.0(2) 


0.2 


0.9(1) 


3.2(1) 


1.4(2) 


2.0(2) 



The density of rungs occupied by two spins at low tem- 
peratures where all holes are bound in HPs is 1 — nh but 
exciting QPs reduces the number of such rungs by one 
for each QP so that the rung density is then 1 — n^ — riQp. 
Our approach to a model for xm is simply to scale the 
form for undoped ladders proposed by Troyer et al. B] 
by this two-spin rung density leading to 



Xm = (l-n h - n QP )(3- 
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M 



where z(f3) — dk{2-K) 1 exp(- 
[A^+4a M (l+cosl)] 1 /2 ||. 

The parameters obtained by a fit of this model to the 
TDMRG data are shown in Tab. || The main change 
upon doping is the decrease of the magnon gap Am 
fll(|[l2|Jl3[ due to interactions between the magnons and 
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TDMRG can be successfully used to reliably simulate 
strongly interacting fermions over a wide temperature 
range. 

Finally, very recently Rommer and Eggert report 
TDMRG calculations for a spin chain with impurities 
|| which uses the same method as we do to overcome 
the problems caused by complex eigenvalues, but they 
did not introduce the re-biorthogonalization which for 
fermionic models, we find is essential. 

We wish to thank E. Heeb, M. Imada, and X. Wang 
for useful discussions. Most of the calculations have been 
performed on the DEC 8400 5/300's of the C4-cluster at 
ETH Zurich. 



FIG. 4. Entropy density of the doped t-J ladder with 
J = t/2 = t'/2 = J'/IO and different hole-densities rih- 



QPs. Due to hybridization with higher lying bands the 
QP bandwidth oqp is also reduced from the leading order 
perturbation result oqp = 2t, but the QP gap Aqp is in 
reasonable agreement with the second order perturbative 
estimate of 0.98L The increase of Aqp (or equivalently 
the binding energy Eb) with nu, can be attributed to an 
effective repulsion between the QPs and HPs. A similar 
increase of the QP gap Aqp was found in Ref. |l2|. This 
is an issue which warrants further investigations. 

Finally, in Fig. |^ we show the entropy density s, 
remapped in the same way to the canonical ensemble 
of fixed n/j. In the limit of T — > oo Sqc = (1 — rih) In 2 — 
rihhinh — (1 — rih) ha(l — rih)- At T/t = 20 the entropy 
has acquired between 99.4% of its maximal value Sqq for 
rih = 0.025 and 99.7% for rih — 0-2. Below the magnon 
gap A a/ , the entropy decreases exponentially, for the un- 
doped Heisenberg ladder down to s = 0. In the pres- 
ence of hole doping, the exponential decrease shows a 
crossover to a linear decrease at low temperatures, as is 
expected for Luther-Emery liquids. This behavior is con- 
sistent with the hard core boson model proposed for the 
magnons. Quantitative fits are however better performed 
on s(fx,T), as we did earlier, due to added uncertainties 
arising from the remapping to constant hole doping. 

In conclusion, we have developed the first numerically 
stable TDMRG algorithm for fermionic systems. This 
has enabled us to calculate accurate results for the mag- 
netic susceptibility x an d the entropy density s of the 
doped t-J ladder with strong exchange on the rungs down 
to low temperatures in the thermodynamic limit of infi- 
nite system size. The system we have studied has two 
crossovers as the spins bind in singlet pairs and the holes 
in hole pairs. These crossovers can be clearly seen in the 
numerical data and demonstrate that this form of the 
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